/**
 * @file prs_max.cpp
 * @author ultralyj (1951578@tongji.edu.cn)
 * @brief 并行max算法的若干实现
 * @version 0.1
 * @date 2021-12-15
 * 
 * @copyright Copyright (c) 2021
 * 
 */
#include "prs_max.h"

/**
 * @brief 求最大值算法的简单实现
 * 
 * @param seqData 单精度浮点数随机序列输入
 * @param seqLength 浮点数序列的长度
 * @return float_t 最大值输出
 */
float_t prs_max(const float_t *seqData, const size_t seqLength)
{
    float_t maxData = -INFINITY;
    for (int64_t i = 0; i < (int64_t)seqLength; i++)
    {
        if (maxData < logf(sqrtf(seqData[i])))
            maxData = logf(sqrtf(seqData[i]));
    }
    return maxData;
}

/**
 * @brief 使用openmp并行算法加速求解最大值
 * 
 * @param seqData 单精度浮点数随机序列输入
 * @param seqLength 浮点数序列的长度
 * @return float_t 最大值输出
 */
float_t prs_max_omp(const float_t *seqData, const size_t seqLength)
{
#define procCnt 16
    /* max变量赋初值 */
    float_t maxData = -INFINITY;
    float_t maxDataThread[procCnt];
    for (size_t i = 0; i < procCnt; i++)
    {
        maxDataThread[i] = -INFINITY;
    }
/* 分16块各自求最值 */
#pragma omp parallel for
    for (int64_t i = 0; i < (int64_t)seqLength; i++)
    {
        int32_t id = omp_get_thread_num();
        if (maxDataThread[id] < logf(sqrtf(seqData[i])))
            maxDataThread[id] = logf(sqrtf(seqData[i]));
    }
    /* 求16个最值中的最值 */
    for (size_t i = 0; i < 16; i++)
    {
        if (maxData < maxDataThread[i])
            maxData = maxDataThread[i];
    }
    return maxData;
}

/**
 * @brief 使用AVX指令集进行并行加速求最值
 * 
 * @param seqData 单精度浮点数随机序列输入
 * @param seqLength 浮点数序列的长度
 * @return float_t 最大值输出
 */
float_t prs_max_avx(const float_t *seqData, const size_t seqLength)
{
#define BlockWidth 8
    float_t maxData = -INFINITY;
    const float_t *procData = seqData;
    __m256 procData_avx = _mm256_setzero_ps();
    float_t *q = (float_t *)&procData_avx;
    /* 8个单精度浮点数同时处理 */
    for (int64_t i = 0; i < (int64_t)seqLength / BlockWidth; i++)
    {
        procData_avx = _mm256_log_ps(_mm256_sqrt_ps(_mm256_load_ps(procData)));
        /* 在8个计算结果之间求最大值 */
        for (int64_t j = 0; j < BlockWidth; j++)
        {
            if (q[j] > maxData)
            {
                maxData = q[j];
            }
        }
        /* 数据向后位移8个float */
        procData += BlockWidth;
    }
    return maxData;
}

/**
 * @brief 使用openmp并行处理分块的最值，并用AVX指令集加速
 * 
 * @param seqData 单精度浮点数随机序列输入
 * @param seqLength 浮点数序列的长度
 * @return float_t 最大值输出
 */
float_t prs_max_omp_avx(const float_t *seqData, const size_t seqLength)
{
#define procCnt 160
    float_t subMax[procCnt];
    float_t max = -INFINITY;
#pragma omp parallel for
    /* 对原序列进行分块处理 */
    for (int32_t i = 0; i < procCnt; i++)
    {
        subMax[i] = prs_max_avx(seqData + seqLength * i / procCnt, seqLength / procCnt);
    }
    for (size_t i = 0; i < procCnt; i++)
    {
        if (max < subMax[i])
        {
            max = subMax[i];
        }
    }
    return max;
}

float_t prs_max_merge(const float_t op1, const float_t op2)
{
    return (op1 > op2) ? op1 : op2;
}